function z = terminalStateCondition(lambda0, xi, ti, tf, m)
% Integrate the system forward in time
[t,s] = ode23(@maxRadiusOrbitTransferEqAndAdjointEq, [ti tf], [xi;lambda0]); 
r_tf = s(end,1);
u_tf = s(end,2);
v_tf = s(end,3);
lambda_r_tf = s(end, 4);
lambda_u_tf = s(end, 5);
lambda_v_tf = s(end, 6);
mu = zeros(3,1);
mu(1) = u_tf; % final state constraint 1
mu(2) = v_tf-1/sqrt(r_tf); % final state constraint 2 
mu(3) = lambda_r_tf+1-lambda_v_tf/(2*r_tf^(3/2)); % final state constraint 3

z = mu - m; % since fsolve solves f(x)=0
